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ABSTRACT 

We show that bars in galaxy models having halos of moderate density and 
a variety of velocity distributions all experience a strong drag from dynamical 
friction unless the halo has large angular momentum in the same sense as the 
disk. The frictional drag decreases the bar pattern speed, driving the corotation 
point out to distances well in excess of those estimated in barred galaxies. The 
halo angular momentum required to avoid strong braking is unrealistically large, 
even when rotation is confined to the inner halo only. We conclude, therefore, 
that bars are able to maintain their observed high pattern speeds only if the 
halo has a central density low enough for the disk to provide most of the central 
attraction in the inner galaxy. We present evidence that this conclusion holds 
for all bright galaxies. 

Subject headings: galaxies: evolution — galaxies: halos — galaxies: 
kinematics and dynamics — Galaxy: halo — Galaxy: structure 

1. Introduction 

The flatness of disk galaxy rotation curves, particularly outside the optical radius, is 
generally interpreted as evidence that they are embedded in massive dark matter (DM) 
halos. Since the mass-to-light ratio of the visible material is uncertain, the central DM 
density is also uncertain; even the best determined 1-D rotation curve is consistent with 
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almost any disk mass up to a maximum that does not over-fit the inner part (van Albada et 
al. 1985). This degeneracy introduces a serious uncertainty into studies of galaxy formation 
and evolution. Here we present an argument that DM halos have the lowest possible central 
density consistent with not being hollow. Following van Albada & Sancisi (1986), we refer 
to such minimum halo models as "maximum disks." 

Strong bars are seen in optical images of roughly 30% of all disk galaxies (Sellwood & 
Wilkinson 1993) and this fraction rises to over 50% in the near IR (Eskridge et al. 2000). 
Their presence makes them ideal probes of the dynamics of the central regions. The rate of 
rotation of a bar can be characterized by the ratio TZ = Dl/cib, where Dl is the corotation 
radius and a# the bar semi-major axis. More precisely, -Dl is the distance from the center 
to the Lagrange point on the bar major axis where the gravitational attraction balances 
the centripetal acceleration in the frame rotating with the bar. Theoretical arguments 
(Contopoulos 1980) require TZ > 1, and there is a prejudice that TZ k, 1 (see Sellwood & 
Wilkinson 1993 for a review). We describe bars for which 1 < TZ < 1.4, i.e. those for which 
corotation is not far beyond the bar's end, as fast. 

The number of barred galaxies with measured 1Z is still quite small (Elmegreen 1996), 
since it requires knowledge of the bar's pattern speed, fl p , which is hard to determine. 
Tremaine & Weinberg (1984a) show that Q p can be measured directly from observations 
of a tracer component that satisfies the continuity equation. To date, their method has 
been applied to just two galaxies: Merrifield & Kuijken (1995, see also Kent 1987) find 
TZ = 1.4 ± 0.3 for NGC 936 while Gerssen et al. (1999) find TZ = 1.15±g;|| for NGC 4596. 
A third case should be completed soon (Debattista & Williams 2000). A less direct, but 
probably reliable, determination of TZ has been made for those few galaxies in which high 
spatial resolution, 2-D gas kinematics have been modeled. Three cases are: TZ = 1.3 in 
NGC 1365 (Lindblad et al. 1996), TZ = 1.3 for NGC 1300 (Lindblad & Kristen 1996), and 
TZ = 1.2 for NGC 4123 (Weiner et al. 2000). A more general argument can be made on 
the basis of the shapes of dust lanes that TZ — 1.2 ± 0.2 (van Albada & Sanders 1982; 
Athanassoula 1992). Other methods to determine TZ rely on identifying resonances, such 
as rings in the disk (e.g. Buta 1986; Buta & Combes 1996), but are less reliable. Thus 
the meager reliable measurements are mostly for galaxies of earlier type, but all indicate 
bars are fast; this conclusion becomes much stronger, and can be extended to later Hubble 
types, if the dust lane argument holds. 

Weinberg (1985) predicted that dynamical friction should brake the rotation rate of 
a bar on a time scale short compared with the ages of galaxies if a substantial density of 
dark matter is present in the region of the bar. In a previous paper (Debattista & Sellwood 
1998, hereafter Paper I), we confirmed this prediction for non-rotating halos with isotropic 
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velocity distributions and concluded that DM halos must have low central densities if 
bars are to remain as fast as those observed. Here we extend this result to halos with 
anisotropic velocity dispersion tensors both with and without net rotation and show that it 
is not significantly altered for any reasonable velocity distribution of the dark matter (cf. 
Tremaine & Ostriker 1999). 



2. Methods 

We present fully self-consistent, 3-D simulations of barred disks embedded in live halos. 
Our simulations start from axisymmetric disk and halo models which are designed to be 
unstable to the formation of bars. 

We adopt the Kuz'min-Toomre disk (Kuz'min 1956; Toomre 1963; Binney & Tremaine 
1987, §2.2.1) because its density drops off less steeply than an exponential of the same scale 
length. We give the disk a Gaussian density profile in the vertical direction with a scale 
height Zd and sharply truncate it at some radius R t . Thus the disk density in cylindrical 
coordinates is 



p(R,z) 



(2n)^z d Rl [1 + (R/R d y]^ R ~ Rt , (i) 



I R > R t 

where R d is the length scale of the disk and f d is the fraction of mass in the disk. 

We set up the initial halo from a distribution function (DF) that is computed to be in 
equilibrium in the presence of the adopted disk, as described in Appendix A. We adopt a 
lowered polytrope form for the DF: 

f(x,v) = CT{E) =c{[-2E(x,v)] n -i - {-2E max ] n ^} (2) 

where C is a normalization constant and n and E max are free parameters. We emphasize 
that the halos generated this way are not polytropes; in particular, the density reaches zero 
at a finite radius for all n > |, whereas true polytropes of index > 5 have non-zero density 
everywhere. By setting E max to the potential energy in the plane of the disk at some radius 
within the grid, we guarantee that no particles are initially on orbits that would take them 
off the grid. When T is a function if E only, the DF is isotropic and the halo almost 
spherical, but in later sections we make T a function of a combination of E and J z which 
leads to anisotropic DFs and spheroidal density distributions. The procedure for selection 
of particles from a DF is described in Appendix B. We set some of our halos into rotation 
by flipping the sign of L z for some fraction of the halo particles, as described in Appendix 
C. In some cases, we flip particles in the inner part of the halo only, with a rule (equation 
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U2 ) that depends on energy in such a way that rotation declines continuously to zero near 



some (spherical) radius r rot . 

We give the disk particles some random velocities with a radial dispersion, a u , set 
to yield a constant Toomre Q, neglecting any corrections arising from disk thickness and 
force softening. We then use the epicyclic approximation to set the azimuthal velocity 
dispersion, a v and the equation for asymmetric drift to set the mean orbital speed. (This 
equation sometimes has no solution at small radii, particularly for large Q, in which case, 
we reduce a u .) Finally, we set the vertical velocity dispersion, a w , from the 1-D vertical 
Jeans equation. This disk set-up procedure is approximate, particularly for larger Q, but in 
practice the disk quickly adjusts to an equilibrium. 

We have generally chosen Q = 0.05 initially in order to hasten the formation of the bar, 
since we are here most interested in the evolution after this event. Models run with higher 
initial Q are not qualitatively different; even though the bars are initially weaker and buckle 
out of the plane more, they end up somewhat longer and are braked to an even greater 
extent. The evolution is not significantly affected by changes to the truncation radius or by 
doubling or halving the initial thickness of the disk. 

We imposed an initial three-fold symmetry on the models by replicating particles in 
sets of six in such a way as to ensure that the total momentum and components of the total 
angular momenta about the x- and y-axes (the z-axis being the symmetry axis) were all 
zero. This prevented the model from rotating or drifting relative to our grid, which could 
lead to excessive and asymmetric loss of particles from the grid in our very long runs. 

We adopt units in which G = M unit = = 1, where G is Newton's constant. The 
total of the disk (M d = f d M tot ) and halo (M h ) masses is M tot = [1 - (1 + iJf/i^)- 1 / 2 ] M unit . 
Our unit of time is therefore (i?d/GM unit ) 1//2 , and our adopted time step is 0.05 in this unit. 

We employ the 3-D Cartesian particle-mesh code described in Sellwood Sz Merritt 
(1994), which uses an efficient FFT-based Poisson solver (James 1977). This choice of 
code does limit us to computing the evolution within a fixed volume, and since we wish to 
retain reasonable spatial resolution within the disk, we have generally been forced to bound 
our halos at a small radius. We find this a reasonable price to pay, since e.g. treecode 
simulations of models with the same number of particles would run ~ 100 times more 
slowly (Sellwood 1997) which would preclude the extensive exploration of parameter space 
we report here. Even the two simulations reported in §|j|, which use larger grids to permit 
more extensive halos, run only ~ 12 times more slowly than our standard grid, and are 
therefore still decisively less expensive than a treecode. These performance comparisons 
are all based on a single processor general purpose workstation; the advantage of a grid 
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Halo type 
Disk type 
Disk scale length, R d 
Disk scale height, z d 
Disk truncation radius, R t 
Disk mass fraction, f d 
Disk Toomre Q 
Halo truncation radius, 

Grid size and type 
Number of particles, N 
Time step 



Lowered polytrope, n = 3 
Kuz'min- Toomre 
5 mesh spaces 
O.LRa 

4# d 
0.3 
0.05 
l2.6R d 

129 x 129 x 129 Cartesian 
300,000 (equal mass) 



0.05 {R 3 JGM U 



1 1/2 



Table 1: Properties and parameters in the canonical simulation. These are default values 
for all simulations reported in this paper, except where noted otherwise. 



code would be even greater on parallel computers where the field evaluation is more easily 
optimized. Athanassoula et al. (1998) give performance comparisons between the grid code 
and a machine with special-purpose hardware (GRAPE3), which clearly must depend upon 
the workstation adopted for comparison, but from their Figure 3, it can be seen that the 
grid code is competitive with their GRAPE3 machine for the grid size and particle number 
we employ. 

The numerical parameters in the simulations we report here are summarized in Table [TJ 
We have checked that our main conclusions are insensitive to changes in particle number, to 
an increase in spatial resolution or the method to determine the gravitational field, etc. In 
particular, the evolution of the pattern speed and bar length on finer grids, or with different 
geometry tracked that on our standard grid pretty well, see Debattista (1998) for details. 
We repeated three simulations with identical physical properties but different random seeds 
in the generation of particles and report the results in Figure [7[ Some particles were lost 
from the grid (typically no more than 3% in the longest runs) almost all of which were halo 
particles. Naturally, these weakly bound particles carried away more than their fair share 
of angular momentum, which decreased by as much as 5% in the worst case. A test on a 
larger grid showed that the evolution is imperceptibly affected by this loss. 



2.1. Pattern speed and Lagrange point 



In order to determine the bar pattern speed, f2 p , and other properties of our models, 
we need a well-defined center about which to perform an expansion. Despite our careful 
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set up, the center of our iV-body system wanders from the center of our grid. Following 
McGlynn (1984), we define the function: 



N 



[1 K 
(xj - x ) 2 + (yj - y ) 2 + (zj - z ) 2 \ (3) 

where (xj, yj, Zj) are the coordinates of the jth particle and (xo, yo, zq) those of an expansion 
center. Minimizing oot with respect to xq, yo and zq gives a centroid for the system. Setting 
k — 1 yields the center of mass, distant particles are weighted more when k > 1 while k < 1 
gives a centroid more sensitive to small scales. We adopt k — |, which removes the dipole 
moment. We determine the centroid for the total system of particles, disk and halo; the 
disk and halo centroids typically differed significantly only during the brief interval when 
the bar buckled, when the position angle of the bar was anyway not well defined. We 
obtain an improved estimate of the centroid position from a single Newton iteration of its 
old value every 20 steps, immediately before each analysis step. We have verified, at a few 
selected times, that further refinement results in a change in the position of the centroid by 
& 0.02i? d . 

We expand the instantaneous distribution of disk particles in logarithmic spirals as 

I N d 

A(m, 7, t) = — exp im (ipj + tan 7 In .R,-) , (4) 
Nd 3=1 

where 7 is the angle between the radius vector and the ridge line of the spiral. Here, (Rj, ipj) 
are cylindrical polar coordinates (with respect to the centroid) of the j-th particle at time 
t. The m = 2, 7 = term of this expansion gives the phase and amplitude of the bar 

A(2,0,t) =Ae 2[v , (5) 

where A is the bar amplitude. We calculate the monotonic angular displacement of the 
bar, 0(t), from (p. We estimate its derivative, Q p (t) = <f>(t), by fitting a straight line to 
25 consecutive data points centered at t, which smoothes out rapid fluctuations in Q p and 
yields an error estimate from the standard error in this linear fit. 

Having determined Q p , we are in a position to calculate D^. We determine the effective 
force along the bar major axis in the disk plane (we average the gravitational force from 
both sides of the center, and use the instantaneous value of Q p ); Dl is the distance from 
the center at which the net force passes through zero. We use the uncertainty in f2 p to 
determine that in directly. This procedure is superior to determining the radius at 
which RQ p intersects an axisymmetric rotation curve; we have found that -Dl is larger by 
more than ~ 5% for strong bars when 1Z ~ 1. 
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It should be noted that the value of D L is affected by the rotation curve shape. As the 
bar slows, the distance from the center to the Lagrange point increases more slowly when 
the rotation curve declines than it would do if the rotation curve were flat. Since, in the 
large majority of our simulations the rotation curve does in fact drop continuously from 
the maximum in the disk, our reported values of both Dl and 1Z are underestimates of the 
values they would have in more realistic models with flat rotation curves (see §5). 

2.2. Bar semi-major axis 

A bar is a straight bisymmetric distortion in the density of a disk; even Fourier 
components of the surface density are therefore ideal for tracing the extent of the bar. Thus, 
for example, Lindblad et al. (1996) and Lindblad & Kristen (1996) used the phase variations 
of the even Fourier components to determine the lengths of the bars in NGC 1365 and 
NGC 1300 respectively. Here we adopt a similar approach for determining the semi-major 
axis, ob, of the bars in our simulations, using only the m = 2 Fourier component. This is 
not always an easy measurement, particularly when the disk possesses m = 2 disturbances 
other than the bar, such as spirals, rings surrounding the bar, and outer oval distortions. 

We divide the disk into radial bins of fixed width 0.16 and determine the amplitude 
and phase of the second sectoral (m = 2) harmonic from the particle positions within each 
annular bin. We estimate the errors o" amp and a p h s using Monte-Carlo measurements of the 
phase and amplitude from synthetic samples of particles drawn from a distribution with 
a given m = 2 amplitude. Fitting these data with a 2-D spline then yielded interpolation 
formulae, giving the uncertainties for the number of particles in each annulus and the 
measured m = 2 amplitude. 

Spirals are the easiest to distinguish, since they generally have a different pattern 
speed from the bar. Thus at a fixed radius, the peak of a spiral's azimuthal position is 
usually different from that of the bar. When the inner part of the spiral lines up with the 
bar, however, there is no way to distinguish between it and the bar; measurements of as 
therefore fluctuate at the beat frequency of the bar and spiral pattern speeds (this same 
beat frequency can be seen in measurements of fl p ). This problem becomes less severe as 
the evolution proceeds, because a rising Q causes the spirals to weaken. 

Oval outer disks are often perpendicular to the bar (in this sense, this outer region can 
be considered as an outer ring of the type Ri in Buta's [1986] classification). The number of 
particles in these ovals is often low, and the error bars associated with their position angle 
correspondingly large enough to confuse the measurement of Ob, especially late in the runs 
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when the surface density just beyond the bar's end has been depleted. To counteract this 
tendency, we did not include radial bins in which the error in position angle is greater than 
~ 80°. An example of a mildly oval outer disk can be seen in Figure 

At later stages in the disk's evolution, the spirals often settle to form a ring around 
the bar, as can be seen in Figure |j. Their position (just outside the bar) and orientation 
(usually aligned with the bar) suggest they are inner rings (Athanassoula & Bosma 1985; 
Buta 1986; Buta & Combes 1996). Since these elliptical rings mostly line up with the bar, 
they constitute an additional m = 2 component locked at the bar's pattern speed which 
further complicates identification of the bar end. When the disk is sub-divided into annuli 
and the m = 2 amplitude plotted as a function of radius, the ring appears as an upwards 
bump in the smooth decrease of amplitude from the bar. We adopt, as one estimate denoted 
obi, the last radial point at which this radially binned m = 2 amplitude did not deviate 
from a linear decrease by more than its standard error, <r amp . 

We obtain a second estimate, as2i from the phase of the m = 2 moment of the disk 
binned as for the clbi measurement. We estimate a# 2 as the radius of the outermost bin for 
which the phase is constant to within the standard error in that radial bin, a p h s - 

We define to be the simple average of clbi an d clb2- m practice, Obi tended to 
underestimate our subjective visual impression of the bar semi-major axis (particularly at 
early times), while as2 tended to overestimate it. We found that the average of these two 
quantities did a very good job of estimating as- We generously define the uncertainty in 
ob, which is not a formal error, as half the difference between clbi and as2- Because of the 
formation of rings and ovals, this uncertainty tends to be largest at late times, and can be 
as large as several disk scale-lengths. 

The uncertainty in Dl is always much less than our generous estimates of the 
uncertainty in Ob, which therefore dominates our quoted uncertainty in 1Z. 



2 The presence of these rings in our simulations has important ramifications for the theory of ring formation 
and interpretation, since rings are often considered to be gas phenomena, but our collisionlcss simulations 
have no gas. 
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2.3. Parameterizing rotation curve decompositions 



One estimator of the relative contributions of the disk and halo to the central attraction 
is the parameter 



where V^disk and K,haio are the circular velocities due to the disk and halo respectively. In 
Paper I, we defined r) exp at -R mj exp, the disk-plane radius at which K,disk would peak for 
an infinite, razor thin exponential disk, as appropriate for the model with the extensive 
halo; we then adopted R m = 1.75R d as the nearest equivalent for the Kuz'min-Toomre 
disks. Here, however, we define r] at the true disk maximum for the Kuz'min-Toomre disk, 
R m = 1.41i? d , since we employ that disk in all the simulations reported here. It should be 
noted that even though the rotation curve evolves as our simulations proceed, our values 
for rj are those of the initial model only. 

Here we explore a wider range of models than in Paper I, with a greater variety of 
rotation curve shapes and find, not surprisingly, that a parameter which depends on the 
ratio of rotation velocities from disk and halo at a single radius does not correlate well with 
the degree of braking we observe in our simulations. Furthermore, halo angular momentum 
changes the evolution of 7Z, so that rj is clearly an inadequate predictor of 1Z even for fixed 
rotation curves. 

A parameter which takes into account both the angular momentum in the halo and the 
full rotation curve might do a better job. We define 



where J 2j h,i is the angular momentum of the i-th halo particle about the symmetry axis, 
r is some arbitrary (spherical) cutoff radius for the summation and J Zt d,i is the angular 
momentum of the i-th disk particle. The quantity r(r ) is the difference between the 
maximum possible and the actual angular momentum content of the halo, expressed as a 
fraction of disk's angular momentum, and is zero for a maximally rotating halo. It can be 
thought of as a measure of the capacity of the inner halo to accept angular momentum; it 
depends both on the halo mass distribution as well as its angular momentum content. 




(6) 



r(r„) 



i,r<r \Jz,h,i\ J z,h,i 
^2i,R<4R d Jz,d,i 



(7) 



3. Massive Halo Models 
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Table 2. All runs in this paper 



Run a 


/d 


Q 


n 


H 




L z ,h pars b r) 


T(3) 


D L /a B (1000) c 


^-lmp 


t(lmp) e 


Comments 




Halo type: F(E) 


*4 


0.3 


0.05 


3 


0.0 


0.0 




2.0 


0.63 


6.6/(2.8 ±0.2) 


1.9 ±0.5 


2000 


Canonical run 


20 


0.3 


0.05 


3 


-0.98 


0.0 


34.97 (1) 


2.0 


1.23 


7.1/(2.0 ±0.2) 


3.5 ± 0.3 


1000 




21 


0.3 


0.05 


3 


0.98 


0.0 


34.97 (1) 


2.0 


0.02 


4.2/(2.5 ±0.3) 


1.7 ±0.2 


1000 




22 


0.3 


0.05 


3 


0.33 


0.0 


0.74 (s) 


2.0 


0.47 


5.6/(3.0 ±0.5) 


1.9 ± 0.3 


1000 




23 


0.3 


0.05 


3 


0.66 


0.0 


5.39 (1) 


2.0 


0.27 


4.9/(3.4 ±0.6) 


1.4 ±0.3 


1000 




24 


0.3 


0.05 


3 


-0.33 


0.0 


0.74 (s) 


2.0 


0.79 


6.9/(2.6 ±0.2) 


2.6 ± 0.2 


1000 




25 


0.3 


0.05 


3 


-0.66 


0.0 


5.39 (1) 


2.0 


0.98 


5.6/(2.2 ±0.2) 


2.6 ± 0.2 


1000 




7 


0.3 


1.0 


3 


0.0 


0.0 




2.0 


0.67 


5.6/(3.0 ±0.4) 


3.3 ± 0.3 


2500 




8 


0.3 


1.5 


3 


0.0 


0.0 




2.0 


0.70 


4.5/(3.0 ±0.2) 


3.2 ± 0.2 


3250 




28 


0.3 


0.05 


3 


0.0 


0.0 




2.0 


0.63 


5.7/(2.3 ±0.2) 


1.8 ± 0.6 


2000 


Twin of Run 4 


33 


0.3 


0.05 


3 


0.98 


0.0 


34.97 (1) 


2.0 


0.02 


3.5/(2.5 ±0.6) 


1.7 ±0.4 


2000 


Twin of Run 21 


34 


0.3 


0.05 


3 


-0.98 


0.0 


34.97 (1) 


2.0 


1.23 


7.1/(2.4 ±0.4) 


3.3 ± 0.7 


2000 


Twin of Run 20 


44 


0.3 


0.05 


3 


0.0 


0.0 




2.0 


0.63 


6.5/(3.9 ± 1.0) 


2.5 ± 0.4 


1500 


z d = 0.2R d 


47 


0.3 


0.05 


3 


0.0 


0.0 




2.0 


0.63 


6.3/(2.6 ±0.2) 


2.4 ±0.1 


1000 


z d = 0.05R d 


Halo type: T{E + 


50 


0.3 


1.0 


3 


0.0 


-0.1 




2.4 


0.89 


6.9/(3.0 ±0.2) 


2.3 ±0.1 


1000 


257 2 x 129; r h = 8R d 


54 


0.3 


1.0 


3 


1.0 


-0.1 


full 


2.4 


0.00 


4.5/(3.4 ±0.6) 


1.4 ±0.3 


1250 


r h = 8i? d 


55 


0.4 


1.0 


3 


0.0 


0.1 




1.4 


0.50 


5.6/(3.7 ±0.7) 


2.0 ±0.2 


1750 


Rt = 8Rd 


56 


0.4 


1.0 


3 


1.0 


0.1 


full 


1.4 


0.00 


4.4/(3.6 ±1.0) 


1.2 ± 0.4 


1500 


Rt = 8Rd 



Halo type: T\E + /3(E) jf] 



123 


0.3 


1.0 


3 


0.0 


-0.2(1 


-^) 2 




2.9 


0.67 


7.1/(3.0 ±0.2) 


2.4 ±0.4 


1250 






124 


0.3 


1.0 


3 


1.0 


-0.2(1 


-z? 


full 


2.9 


0.00 


3.5/(2.6 ±0.4) 


1.3 ± 0.1 


1250 






137 


0.3 


1.0 


3 


0.04 


-0.2(1 


-t? 


2.0 


2.9 


0.55 


6.7/(3.3 ±0.2) 


2.1 ± 0.2 


1250 


r ro t 


= 2R d 


138 


0.3 


1.0 


3 


0.15 


-0.2(1 




1.5 


2.9 


0.28 


4.9/(2.8 ±0.5) 


1.7 ±0.3 


1000 


Trot 


= 4R d 


139 


0.3 


1.0 


3 


0.50 


-0.2(1 




1.2 


2.9 


0.06 


4.2/(2.4 ±0.2) 


1.6 ±0.3 


1250 


Trot 


= 6Rd 


142 


0.3 


1.0 


3 


0.0 


-0.5(1 






4.8 


0.67 


5.8/(3.0 ±0.2) 


2.1 ± 0.1 


2000 






143 


0.3 


1.0 


3 


0.11 


-0.5(1 




1.5 


4.8 


0.33 


3.3/(3.0 ±0.8) 


1.5 ± 0.3 


2000 


r ro t 


= 4R d 


144 


0.3 


1.0 


3 


1.0 


-0.5(1 


-t? 


full 


4.8 


0.00 


3.6/(2.8 ±0.3) 


1.3 ±0.2 


2000 






145 


0.3 


1.0 


3 


0.03 


-0.5(1 


-z? 


2.0 


4.8 


0.58 


3.8/(2.2 ±0.2) 


1.6 ±0.2 


2000 


rrot 


= 3R d 


146 


0.3 


1.0 


3 


0.47 


-0.5(1 


-^) 2 


1.2 


4.8 


0.07 


3.5/(2.9 ±0.4) 


1.3 ±0.1 


2000 


rrot 


= 7R d 


147 


0.3 


1.0 


3 


0.04 


-0.5(1 


-c) 2 


1.8 


4.8 


0.51 


3.6/(2.6 ±0.2) 


2.3 ±0.2 


3500 


rrot 
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Fig. 1. — (a) The initial rotation curve of the canonical run. The unbroken line is the total 
rotation curve, the dot-dashed line is the disk contribution and the dashed line is the halo 
contribution, (b) The rotation curve of the axisymmetric (azimuthally shuffled) particle 
distribution towards the end of the simulation. The extra (dotted) line shows the halo 
contribution at t — for comparison. Note that angular momentum redistribution resulted 
in a more concentrated disk, and a slightly lower density halo. 

We begin by describing a set of experiments with disks embedded in moderately dense 
halos. Halos of significantly greater density than we use here would inhibit the formation of 
the bar (Ostriker & Peebles 1973; Toomre 1981). We already demonstrated in Paper I that 
such a halo having an isotropic velocity distribution would brake the bar to an unacceptable 
extent. Here we determine the extent to which braking is affected by giving the halo net 
angular momentum, or an anisotropic distribution of velocities, or both. 



3.1. Canonical simulation 

Our canonical simulation (run 4) is the most halo-dominated model reported in Paper 
I, which we describe in more detail here. The initial rotation curve for this sub-maximal 
disk model (Figure |I]a) drops unrealistically beyond the disk edge because the halo density 
drops smoothly to zero at the grid edge. The initial properties, numerical parameters and 
principal result are given in Tables [I] & 0. 

The disk quickly forms a bar, as shown in Figure @(a). Shortly after its formation, at 
t ~ 150, the bar buckled (Combes & Sanders 1981; Raha et al. 1991) very mildly, causing it 
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3.0/(2.6 ±0.2) 
3.8/(2.4 ±0.2) 


1.6 ± 0.3 2000 Maximum disk run 
2.0 ± 0.4 2000 Control run 



a R,uns also presented in Paper I are marked with an asterisk. 

b Parameter settings in generating halo angular momentum. For isotropic halos, we report the value of a followed by (1) or (s) 



depending upon whether the angular momentum needed is large or small in equation CI, For varying anisotropy halos, we report the 
value of d. An entry of "full" in the column means that all retrograde particles have been flipped to give the initial setup. 

c The value of D^/ag by t = 1000, the minimum duration of all simulations. 

d The value of 1Z at the last measured point, which is not always when the bar has finished slowing down. All bars in our simulations 
form with TZ ~ 1. 

e The time at which 7£i mp is reported. Note that one rotation at i? d takes 19 time units in the canonical run. 



13 




Fig. 2. — The evolution of the canonical run. (a) The amplitude of the bar grows very 
rapidly, experiences a weak buckling around t = 200 and then continues to grow slowly, 
(b) The total (topmost), disk (middle) and halo (bottom) angular momentum. Angular 
momentum is lost from the disk and gained by the halo, with a little carried off the grid 
by escaping particles, (c) The pattern speed drops rapidly soon after the bar forms, but it 
reaches a constant value by t ~ 1600. (d) as (crosses) and Dl (circles). The horizontal line 
joining the last 6 values of Ob shows the weighted average final value of a# , which gives 
K = 2.1 ±0.2. 

to thicken. The continuing slow rise in A after this time results from gradual trapping of 
additional disk particles into the bar, often associated with spiral activity (Sellwood 1981); 
we describe this process as secondary bar growth. Spiral activity gradually declines as the 
outer disk heats to Q ^ 4. Figure |3| shows contour plots of the disk at several instances 
which clearly reveal the thickening of the disk, and the peanut shape of the bar. Figure ^ 
shows the distribution of disk particles at t — 1000; note the prominent inner ring and the 
oval outer disk. 

Once the bar forms, Q p has a well-defined value (Figure ^c) which drops rapidly at 
first, but reaches a constant level towards the simulation's end at t — 2000. Figure ||(d) 
shows that ag increases only mildly, whereas Dl increases rapidly at first, later reaching a 
constant value, reflecting the behavior of Q p in Figure 0(c). 
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Fig. 3. — Contours of projected disk density in the canonical simulation at four instants. 
The bar has been rotated into the x-axis. Contours are logarithmically spaced. The circle 
in the (x,y) projection shows a B = 2.48, 2.80, 3.76 & 4.32 at t = 500, 1000, 1500 & 2000 
respectively. Note that the outer disk is distinctly oval, and the peanut-shaped isophotes in 
the edge-on view. 
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Fig. 4. — The disk particles in the canonical simulation at t — 1000. Note the prominent 
ring just outside the bar and the oval outer disk. The box marks the full extent of the grid. 

Figure |2](b) shows that the drop in Q p is associated with a substantial transfer of 
angular momentum from the disk to the halo. The torque which produces this angular 
momentum exchange arises from dynamical friction on the bar as it moves through the 
halo. The bar induces an m = 2 response in the halo which develops very soon after the bar 
forms (Figure |5]a). The response lags the position angle of the bar (Figure [|b) at first, but 
gradually shifts into alignment with the bar as the torque drops. 

It is interesting that the bar survived such strong braking (cf. Kormendy 1979), 
which reduced its pattern speed by a factor of ~ 5. Most theoretical work, starting 
from Contopoulos (1980, see Sellwood & Wilkinson 1993 for a review) has suggested that 
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Fig. 5. — (a) The amplitude of an m = 2 coefficient of the halo response in the canonical 
run. (b) The phase difference between the bar and the halo response. 

self-consistent bars should nearly fill their corotation circles; our simulation is a clear 
counter-example. Its pronounced butterfly-shape when seen from above may be consistent 
with the prediction by Teuben & Sanders (1985) that slow bars require a large fraction of 
"box" orbits. 

The rearrangement of angular momentum altered the density distributions in both 
the disk and halo, causing the rotation curve to change to that shown in Figure 0(b). 
The central density of the disk rose significantly but the density profile of the halo barely 
changed, despite all the work done on it by the bar (Figure ||b). Such a small change 
underscores how difficult it is for any dynamical interaction with the disk to modify the 
halo density profile. 

We emphasize that a more realistic flat rotation curve model would require a more 
massive and extended halo. Not only might this increase dynamical friction, and slow the 
bar still more, but also the Lagrange point would lie further out in the disk, increasing the 
value of 71; this effect alone would increase D L by more than 30% in this run at t = 2000. 

3.2. Halo rotation 

All our models reported in Paper I, including our canonical model, have isotropic halos 
with no initial net angular momentum. Here we turn our attention to the effects of halos 
with net rotation. At first, we simply change the sign of J z for some halo particles according 
to the rule described in Appendix 0. 
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Fig. 6. — The mean rotation rates (full-drawn lines) of the halo particles in the simulations 
discussed in §|3~2|. The circular speed in the mid plane is shown by the dotted curve. 



We ran a series of experiments, summarized in Table 0, in which the halo angular 
momentum, H = L z ^/ L ZtUiax , was ±0.33, ±0.67 and ±0.98, where L z ^ max is the maximum 
achievable if the angular momentum of every particle is made positive. (The corresponding 
values of the dimensionless spin parameter, A = LIE] 1 / 2 /(GM 5 / 2 ) are 0.05, 0.11 and 0.16, 
respectively). The mean rotation speeds of the halo particles are shown in Figure |6|. Note 
that we did not continue all these simulations until the bar had finished slowing down. 

Figure ^ shows that the value of 1Z reached by t = 1000 correlates strongly with the 
angular momentum content of the halo. As found previously by Athanassoula (1996), the 
bars in models with retrograde halos are more strongly braked, while those in prograde 
halos less so, in comparison with the non-rotating case. Nevertheless, even when direct 
rotation in the halo is maximized, TZ = 1.7 ± 0.4 by the end of the run despite the fact that 
bar was weaker. The value of A in the maximally rotating models settled at some 70% 
of its value in the non-rotating simulation, suggesting that secondary bar growth, with a 
concomitant increase in friction, may be enhanced by angular momentum loss to the halo. 
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Fig. 7. — The variation of Q p and TZ with L z ^/L Ztmax at t = 1000 for all the simulations 
with isotropic halos that were set into rotation. Note that there are two runs at each of 
L Zjh /L 2 max = —0.98, 0.0 and 0.98 giving some measure of the inherent scatter. 

3.3. Anisotropic halos 

We next present simulations with halos having somewhat more general DFs that yield 
anisotropic velocity distributions even in the absence of net rotation. Inspired by Osipkov 
(1979) and Merritt (1985), we chose the form T[E + /3J^), where T has the same form as 
in equation (Q). The halos are oblate and azimuthally biased when (3 < and prolate and 
radially biased when (5 > 0. Anisotropy changes the distribution of mass within the halo; 
we therefore adjusted the halo mass fraction, / h , and/or the halo truncation radius, r ha i , 
to generate models with values of i] similar to that in our canonical simulation. 

Some initial tests revealed that too pronounced an azimuthal bias {(3 < —0.2) led to 
strongly lop-sided halos which interfered considerably with the development of the bar. 
These m = 1 instabilities, which appear to be of the kind discussed by Sellwood & Valluri 
(1997) for counter- rotating oblate spheroids, produced much larger asymmetries than those 
generally observed in real galaxies. We therefore report only those simulations which did 
not become strongly lop-sided. Disk-halo interactions via such m = 1 modes are interesting 
in their own right and deserve a separate study. 

The possible parameter space when the halos are allowed to be anisotropic is very 
large; we considered only two main models, an azimuthally biased model near the limit of 
m = 1 stability, and a radially biased model (Table Q). Their rotation curves are shown in 




Fig. 8. — Rotation curve decomposition of the models having (a) the azimuthally 
biased/oblate halo and (b) the radially biased/prolate halo. The solid line is the total 
rotation curve, the dot-dashed line is the disk contribution and the dashed line is the halo 
contribution. The mean rotation speeds of the halo particles in the fully-rotating models are 
shown by the dotted lines. 

Figure [5[ In both cases, we also tested fully-rotating versions of these models, which should 
have the weakest friction (c/. §|3.2|). 

The bias towards large angular momenta introduced by setting j3 < tends to place 
more halo material at large radii at the expense of small radii and the inner rotation curve 
becomes dominated by the disk. To compensate, we further decreased the radial extent of 
the azimuthally biased halos. The oblate halo has an axis ratio ~ 0.7 at 6 R^- The bar that 
formed in a non-rotating halo was strongly braked, but friction is greatly reduced in the 
run with maximum halo rotation (A = 0.23). 

The halos generated by radially biased DFs have larger densities at smaller radii, all 
other things being equal, than do isotropic halos. We therefore reduced the halo mass 
fraction in order to avoid a system which was too halo dominated. The halo was prolate at 
large radii, having an axis ratio ~ 1.16 at 4Ra, but it becomes oblate at Rd because of the 
disk's influence. Since a trial run showed that secondary bar growth can be quite rapid in 
this model, we extended the disk to R t = 8 R^- Once again, the bar is strongly braked in 
a simulation with no net halo rotation, but remains fast when the halo rotates maximally 
(A = 0.13). 
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Fig. 9. — (a) The halo flattening for (3q = —0.2 (full-drawn line) and for /3q = —0.5 (dotted 
line). These halos vary from nearly spherical at large radii to highly flattened in the inner 
parts, (b) The initial mean orbital speeds of the halo particles (full-drawn curves) in four 
simulations of the @ = —0.2 case with rotating halos (see text). The dotted curve shows 
the circular orbital speed in the disk. 

3.4. Radially varying anisotropy and rotation 

Tremaine & Ostriker (1999) suggest that the stringent limit on the halo central density 
that we obtained for isotropic halos (Paper I) could be relaxed if the halo had significant 
rotation in its inner parts only. They proposed that the inner halo had itself been torqued 
up and flattened by interactions with the disk. 

To test their hypothesis, we create halos with varying anisotropy. We use polytrope-like 
distribution functions as before, T{E + (3 J%), but now we let (3 itself be a function of energy, 

P(E)=p (l-e 2 ) (8) 

where e — (E — E m i n ) / (-E max — E m i n ). Here, E max and E min are the values of the potential 
at the grid edge in the disk plane and the center of the system respectively. We set the free 
parameter (3q = —0.2 in order to generate azimuthally biased, oblate halos, in line with the 
prediction of Tremaine & Ostriker. Figure ^| shows the axis ratio of the halo as a function 



r(3) 

Fig. 10. — The variation of 1Z with T(3) for the (3q = —0.2 cases. The largest value of T is 
for a non-rotating halo and the dotted line marks 7Z = 1.4. The values are for t = 1000 or 
t = 1250 - the later time in simulations in which friction is still acting. 

of radius, which varies from ~ 0.5 at the center to spherical at the edge. 

We ran a set of three experiments (runs 137-139 listed in Table ||) in which we 
introduced rotation by flipping retrograde halo particles with a probability that was 
a function of energy given by equation ( |C2|) , which gives something close to maximal 
rotation in the halo to some limiting energy; beyond some spherical radius r rot the halo is 
non- rotating (Figure |9|b). We varied this critical energy in this set of experiments. Two 
other experiments (runs 123 & 124) bracket them in terms of angular momentum content 
by having, respectively, no halo angular momentum and the maximum possible. 

Again, the results are given in Table |2|. As the angular momentum content of the halo 
is increased, the bar slows down less, but even when r rot = 6 (Figure |^b) the final value of 
7Z is still quite large. The parameter 1Z remains acceptably small only for the fully rotating 
case. 

Mildly rotating halos in a similar set of experiments with a still more flattened halo 
(/?o = —0.5), produced much weaker bars at first, because the bar buckled more violently; 
it seems that a different kind of coupling to the halo occurred in these cases. As these 
weak bars were slowed to a lesser extent after a fixed amount of evolution, we continued 
one simulation to t = 3500, and found that the bar recovered and substantial friction again 
developed leading to a slow bar (7Z = 2.26 ± 0.20). 

The mean orbital speed of the halo particles (Figure ||b) already indicates that 
significant halo rotation out to quite a large radius is needed to avoid strong friction. Figure 
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Fig. 11. — Contributions to the rotation curves from the n = 1.6 halo (dot-dashed), n = 4 
halo (dashed), and the disk only (unbroken line). The curves converge at large radii because 
these distributions all have equal total mass. 



10] underscores just how much halo angular momentum this implies: we see that fast bars 
require T(3) ^ 0.15 which is ~ 0.5 less than the value for the non-rotating halo. Thus the 
halo angular momentum inside r = 3 has to be fully 50% of that of the entire disk. 

Tremaine & Ostriker argue that strong halo rotation could be induced out to ~ 3 kpc 
in a Hubble time, but we have shown here that the halo angular momentum requirements to 
avoid strong friction are considerably greater than their mechanism seems able to produce. 



4. Halo Density and Concentration 

In this section we report experiments in which the halo mass and concentration 
are varied, still with the halos confined to small radii, as in §3. In §H, we discuss more 
realistic models with extended halos. We have varied both the halo mass fraction and the 
polytrope index, n. Increasing n leads to more concentrated halos, resulting in a larger 
halo contribution to the inner rotation curve (Figure |ll|). Table ^| lists the parameters of 
these simulations and gives our principal result. Note that these simulations represent three 
series, with varying halo mass density at fixed n. 
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Fig. 12. — The variation of 1Z with rj for different values of the polytrope index: triangles 
n = 4, squares n = 3, and circles n = 1.6. We plot averages of three values of 7Z near the 
end of the simulations when friction is low, but not always zero. 

The evolution of the n = 3 runs has already been presented in Paper I. We find that 
decreasing the halo contribution to the rotation curve leads to a marked decrease in 7Z. 
A fast bar (with a pattern speed that continued to decrease very slowly) survived in the 
simulation with the least massive halo of this series. 



Figure |T| summarizes the values of 1Z obtained from models with different polytrope 
indices. Although there is no trend when all the points are taken together, within each 
series of runs at fixed n, bars with larger r\ end up faster. Note that it would be incorrect 
to conclude from this Figure that less braking occurs for a given r\ as the central density 
increases. This apparent trend results from two different effects: First, increasing n for 
fixed halo mass leads to greater halo concentration, depleting halo material at larger radii 
in our halos (which we confined to a small volume), thereby reducing friction somewhat. At 
the same time, a more sharply peaked halo rotation curve leads to a smaller rj at fixed Mh, 
and also makes the rotation curve drop more steeply, reducing D L . 
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Fig. 13. — The unbroken curves show the circular speeds of the maximum disk (bold) and 
the control experiments. The dot-dashed line shows the disk contribution (common to both 
experiments) and the dashed lines are the halo contributions. Note that both rotation curves 
are almost flat out to 



We next describe two models in which the halo extended to r = 25 Rd, which is large 
enough to achieve a flat rotation curve with moderate to low central densities in the halo. 
We had to increase the grid to 25 7 3 - an eight-fold increase over that used in most of the 
above experiments, making these runs much more expensive. A large parameter space 
search, such as that described in the previous two sections, would be prohibitively expensive 
with this larger grid. 

We already reported a maximum disk model with an extensive halo in Paper I whose 
evolution was computed using a cylindrical polar grid. One of the two models presented in 
this section closely resembles it, but has a different disk type and is run on a Cartesian grid. 
The other model discussed here is a "control" experiment with a similar rotation curve but 
with a somewhat denser halo. 

The disk, which was truncated at 8 Rd, accounts for 17% of the total mass. The 
polytrope index n = 2 for the maximum disk model, whereas n = 3 for the control run, this 
being the only difference between the two simulations. The resulting rotation curves, shown 
in Figure [13], are quite flat out to the disk truncation radius. Both have a substantial disk 
contribution at small radii: r\ = 7.0 in the maximum disk case and i] = 3.8 in the control 
experiment. 



5. 



Models with more Extensive Halos 
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Fig. 14. — The upper panels show the evolution of f2 p for the maximum disk and the control 
simulations while Dl (circles) and a# (crosses) are shown the lower panels. The maximum 
disk is shown on the left (a & c) and the control run on the right (b & d). 

The bars grew very rapidly and did not buckle much, reaching an amplitude some 
~ 5% lower than in our canonical run. Most of the angular momentum lost by the inner 
disk after the bar forms ends up in the halo for both runs, but the outer disk continues to 
accept some of the bar's angular momentum. 



Figure [14] shows the evolution of fi p , a# and Dl for both simulations. The bar in the 



control simulation is slow already by t ~ 1500 and 1Z = 1.98 ± 0.35 by the end. The bar in 
the maximum disk simulation, on the other hand, is acceptably fast at the end of the run, 
with 1Z = 1.57 ± 0.27, but only barely so and it is continuing to slow. We have therefore 
identified a region of parameter space where fast bars can survive for more than ~ 30 
orbital periods at R = 1.4 (just outside the half-mass radius of the disk). 
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6. Synthesis 

We now seek some way to synthesize all these numerical results. We first show that 
the frictional torque acting on the bar from the halo behaves in some sense as might be 
expected from standard dynamical friction. However, the total angular momentum loss 
which occurs before the halo response locks into phase with the bar can be determined only 
numerically; we attempt to relate the final pattern speed of the bar to the ability of the 
halo to accept angular momentum. 



6.1. The Chandrasekhar formula 

Chandrasekhar's (1943) demonstration of a secular drag force on a massive perturber 
moving in a straight line through a uniform, infinite background sea of low mass particles 
differs in many significant ways from the rotation of a bar through an inhomogeneous halo. 
It is now part of received wisdom that his formula works better for a perturber moving 
through an inhomogeneous system, using local values of the density etc., than we have any 
right to expect (Binney & Tremaine 1987, §7.1). We here show that the frictional torque 
also behaves very roughly in accordance with his formula for a rotating bar - at least over 
the period after the bar has formed and settled and before the halo response becomes 
aligned with the bar. 

In the limit when the perturber's mass M is much larger than the masses of the 
background particles, Chandrasekhar's formula for the frictional force is 

,dv A An In AG 2 pM 2 fv\ . . 

M— = -v g(-), (9) 



dt a 2 \o\ 

where v is the velocity of the perturber, v is a unit vector in the same direction and p and 
cr are respectively the density and velocity dispersion of the background. The v 2 factor in 
the denominator of formula (7-18) of Binney & Tremaine has here been replaced by a 2 in 
order to subsume all the velocity dependent factors into the dimensionless function g, which 



is shown in the Figure 15(a) for the case when the velocity distribution of the background 
particles is Gaussian. It can be seen that friction increases as the speed of the perturber 
rises, reaching a maximum when v ~ 1.37er and then decays monotonically for all higher 
speeds. 

The rate of gain of angular momentum of the halo in our simulations is clearly the 
frictional torque t z on it arising from the bar. Note that this measurement is independent of 
bar's reaction to the loss of angular momentum, and therefore does not involve, for example, 
any estimate of its effective moment of inertia. The halo torque could be interpreted as the 




Fig. 15. — (a) The function g denned in equation @. (b) The scaled halo torque T z 
(equation pi]) measured from many simulations all having the same disk/halo mass. See text 
for an explanation of the line. The symbols are as follows: canonical model (squares), hot 
disk models (circles), retrograde halos (triangles), prograde halos (pluses), other polytrope 
indices (crosses), rotating halos with radially varying anisotropy (diamonds). 

frictional force times some characteristic lever arm of length R. We therefore plot in Figure 



151(b) the quantity 



(r> 2 

z R(A) 2 p 1 Uj 

against the speed of the bar perturbation through the background halo at that radius, 
v' = RQ P — (v), normalized by the halo velocity dispersion. The running averages in this 
expression are over 50 time units, and we include A to take account of the variations in bar 
mass in these equal mass disks. The scaling of the ordinate is arbitrary, therefore. We adopt 
R = 3i? m= 2,max (i-e. three times the radius where the m = 2 Fourier component peaks); 
other values of R show the same general trend, but we found that this choice minimized the 
scatter. We adopt local values for p and a (averaged over the range < R < 3i? m=2i max), 
and we use cr 2 = a\ + cr? + a 2 z . 



We plot a curve from one simulation in Figure |15|(b) and a number of points from 
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other simulations. The curve shows the entire evolution of T z , from the moment the bar 
first forms to the end when friction is very low, for a maximally retrograde model (run 20). 
The time evolution along this line is from right to left and can be described as follows: (1) 
the initial spike occurs as the halo response develops soon after the bar forms, (2) the curve 
then dips as the bar buckles, (3) after which the line follows roughly the trend indicated by 
the points as the bar is braked steadily, and (4) it drops down to low values as the halo 
response locks into alignment with the bar. The fluctuations in T z in this run are fairly 
typical; they are larger in some models and less in others. 

We should not expect Chandrasekhar-style friction in any part of this evolution except 
for period (3) after a quasi-steady halo response is established, and before orbit trapping 
becomes significant. Thus we have tried to include points in this plot from other runs 
during the steady friction period, although we obviously failed for the cluster of points near 
v'/a = 0.6 and T z > 0. 

Although there is considerable scatter in these measurements, the general rise over the 
range of the abscissae is similar to the theoretical curve in Figure |I~5|(a), although our data 
show only a rising trend. We are encouraged that a trend shows through despite the crude 
approximations we adopt: we identify a single radius whereas the entire halo probably 
contributes, not all our bars have identically the same density profile, our polytrope halos 
do not have a precisely Gaussian velocity distribution, etc. While perhaps not entirely 
convincing, this Figure does suggest (i) that the torque is indeed from the usual dynamical 
friction, (ii) that most of the drag seems to arise from the region just beyond the end of the 
bar, and (iii) that the characteristic speed is generally less than the halo dispersion. 

It should be noted that the velocity dependence in Figure 0(b) is the opposite of that 
predicted by Weinberg (1985), who suggested that friction would increase if the halo was 
made to rotate in a prograde sense; one interpretation of Weinberg's prediction is that in 
his case most of the friction arises from a perturber speed v ^ lAa. The fact that we find 
the opposite behavior may result from his assumption of an infinite isothermal sphere for 
the halo, whereas our halos have a very limited radial extent. 

6.2. Constraints from T 

In our simulations with strongly confined halos, friction seems to drop to zero before 
the bar is brought to rest relative to the streaming speed of the halo particles. The torque 
vanishes when the induced bi-symmetric distortion in the halo locks into alignment with 
the bar. This locking effect appears to be the result of non-linear trapping of orbits, a 
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Fig. 16. — The distribution of 43 simulations in the T(3), 1Z plane. 



process described as "dynamical feedback" by Tremaine & Weinberg (1984b). Note that 
this locking phenomenon did not occur in our simulations with more extended halos, where 
braking persisted for as long as we ran them (see Figure IT4T). 



We have searched for a predictor of the final bar pattern speed, but have not found 
anything simple - perhaps because none exists. The best we have come up with is the 
parameter Y introduced in § [2.3| . After some experimentation, we found T(3), evaluated 
using the initial values of J z ^ and J z ^, to be the most useful. Figure [16| summarizes the 
last measured values of TZ from an ensemble of 43 different simulations plotted against 
T(3). These simulations include various different polytrope indices, different anisotropies, 
different halo rotations (and different distributions of halo angular momentum), different 
Toomre Qs, different disk thicknesses and different halo masses. Not all of these simulations 
have been evolved to a steady state; we generally stopped the simulation once we found the 
bar to be slow (which may account for some of the scatter in the distributions). 

We draw the following conclusions from the rising trend in this Figure: (1) Bars can 
generally remain fast when T(3) < 0.4, although some are significantly braked. (2) When 
T(3) ^ 0.4, bars generally end up slow. Strong braking can be avoided when T(3) 0.4 
when either the model has a rapidly dropping rotation curve, or a flattened halo with most 
of its mass outside the bar region.^ 

Our parameter T is the best we have been able to find to predict the value of TZ. A 



3 The positions of models with radially varying anisotropy in this plot is more than usually sensitive to 
the value of ro used in computing T. 
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range of values of T can be determined for any rotation curve fit, but, being dark, no value 
of r can be pinned on a halo. Constraints on the actual angular momentum content of dark 
halos might be possible from studies of the faint luminous halos that may trace the rotation 
of the dark halo (being subject to similar dynamics). The Tremaine-Ostriker hypothesis 
calls for low T due to high rotation in the halo. We favor low T due to low halo mass 
fraction (we have argued that the Tremaine-Ostriker hypothesis may need rather unlikely 
levels of angular momentum in the inner parts of halos). The reality may be somewhere in 
between these two limits. 

7. Discussion 

7.1. Comparison With Real Bars: NGC 936 

We need to show that our iV-body bars are similar in strength to bars in real galaxies. 
A photometric comparison would not be conclusive because the strongly non-axisymmetric 
light distribution in a galaxy may not reflect the true distribution of mass. We have 
therefore attempted a kinematic comparison using data from NGC 936, a well-studied 
SBO galaxy having sufficient published data on the stellar velocity field for our purposes. 
Its other advantages are that it appears to be relatively dust-free and is known to have a 
fast bar (Kent 1987; Merrifield & Kuijken 1995). The galaxy has a dense bulge (Kent & 
Glaudell 1989), however, unlike in our simulations. 

We adopt Kormendy's (1983) estimates of the projection geometry, the bar position 
angle, the de-projected bar semi-major axis and use his velocity measurements at five 
different slit position angles. We scale our models by setting Ob = 50" and rotate and 
project them as we view NGC 936. We then compute the mean projected line-of-sight 
velocity (Vi os ) of the particles in the model, averaging the approaching and leading sides to 
maximize the number of particles in our pseudo-slits, and determine the velocity scaling by 
minimizing x 2 between the observations and the model. When making this comparison, 
we use data in the circular annulus (in the galaxy plane) 0.6 < Rjas < 1-2 to exclude the 
region dominated by the bulge and the disk well outside the bar. 

Our canonical model at early times compares very well with NGC 936. We find reduced 
X 2 = 0.70 (for 22 degrees of freedom) at t = 250, and a visual comparison of V[ os shows that 
the variations with position angle track those in the galaxy very well, as they must do for 
this good a fit. (For comparison, we obtain a reduced x 2 = 1-96 if we erase the bar from 
our model by randomizing the azimuthal positions of the particles.) At later times, after 
the bar has been slowed by a factor ~ 5, we find reduced x 2 — 5.94, again showing that 
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NGC 936 is grossly inconsistent with a slow bar. 

The maximum disk model is not quite as impressively similar to NGC 936: reduced 
X 2 = 1-20 at t — 350 which is still acceptable, but reduced \ 2 = 1.54 at t — 1150, which is 
marginally so. The fits to "control" run are again worse: reduced \ 2 = 1-61 at t — 318, 
when the bar is still fast, and reduced % 2 = 1.75 at t = 2000 when the bar has slowed. In 
both cases, the absence of a massive central spheroidal bulge may be partially responsible 
for the poorer fits. 

Thus the bars in our simulations are quite similar to that in NGC 936 when they first 
form, but are clearly inconsistent when they have been strongly braked. 

7.2. Neglected Physics 

Our simulations are of the stellar and DM components of a barred galaxy and do not 
include any gas. It is well known that gas flows in barred potentials produce large-scale 
shocks offset from the bar major axis (e.g. Athanassoula 1992). The asymmetric gas 
distribution loses angular momentum to the bar and gas flows towards the center. The 
small gas fraction in most galaxies, together with the relatively short lever arm, mean that 
the angular momentum given to the bar by gas could not possibly compete with that lost 
to the halo through dynamical friction - e.g. friction removed ~ 40% of the disk's angular 
momentum in our canonical simulation (Figure 0). Furthermore, gas-poor SB0 galaxies, 
such as NGC 936 and NGC 4596, have fast bars. 

Gas inflow has a second effect, however: it deepens the gravitational potential at the 
bar center causing the bar to speed up slightly. This can be a small effect at most, since 
excessive inflow will destroy the bar. The total mass accumulated in the center cannot 
exceed 5% of the disk mass (Norman et al. 1996), and is probably much less; even in this 
extreme case, the increase in bar pattern speed was only some 40% (Sellwood & Debattista 
1996, Figure 6). 

A bulge component inside the inner Lindblad resonance might act as a source of angular 
momentum for the bar. Bulges can be quite massive (e.g. NGC 936, Kent Sz Glaudell 1989) 
and are often in rapid rotation (e.g. Kormendy & Illingworth 1982), but their small radial 
extent limits their angular momentum content. It is possible to think of the inner parts 
of the rapidly rotating halo in some of our simulations as representing a bulge. The bar is 
still strongly braked in such cases (§|3.2|), suggesting that not even rapidly rotating bulges 
can prevent bar slow-down in sub- maximal disks. Weinberg (1985) suggested that the inner 
disk could also be a source of angular momentum for the bar, but we have found that bars 
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always lose angular momentum to the disk, not the other way round. 

We have not included the effects of late gas infall onto the disks which would enhance 
secondary bar growth (e.g. Sellwood & Moore 1999). Irrespective of the rate of bar growth 
by this mechanism, it inevitably leads to increased friction, making it still more difficult for 
the value of 1Z to decrease. 

In addition to massive gas inflow, bars can be destroyed by satellite impacts. Our 
simulations do not take either process into account. If these processes are to account for 
the absence of slow bars, they would have to act efficiently and quickly, and new bars 
would have to form again to maintain the observed high fraction of barred galaxies. The 
rapid formation of a new bar is difficult to arrange, since bar destruction leaves the disk 
dynamically hot, and the inflow destruction mechanism gives rise to a more steeply rising 
inner rotation curve. Both factors limit the disk's ability to form new bars, by making it 
less responsive and by cutting the feed-back loop to the swing amplifier (Toomre 1981). The 
revival of a bar after one has dissolved requires the accretion of so much dynamically cool 
material (Sellwood & Moore 1999) that it is unlikely to occur more than once in a galaxy. 

While the bars in our simulations are comparably strong to that in a real galaxy 
(§ [7.1|) , a systematic difference with early-type galaxies is the absence of a dense bulge in 
our simulations. It is possible that simulations with dense bulges behave differently, but 
it seems unlikely that they would. Orbit studies (Athanassoula 1992) and simulations 
(Sellwood 1989; Sellwood & Moore 1999) reveal that bars in galaxies with dense centers are 
supported by the same orbit family and behave similarly to those in which the center is 
more uniform. 



7.3. Scaling to NGC 3198 

Before discussing the implications of our results for real galaxies, we need to determine 
how they should be scaled. The de facto standard galaxy in the dark matter halo debate is 
NGC 3198; even though it is not a strongly barred galaxy, we nevertheless scale our models 
to the data of van Albada et al. (1985) for this system. 

We use our two extensive halo systems from §|5] for this comparison. Scaling to the 
observed rotation curve determines the length and the velocity scales. We first adopted 
Rd = 3.0 kpc for both simulations, which differs from the exponential scale of 2.68 kpc that 
fits NGC 3198 well (van Albada et al. 1985) since our models have Kuz'min- Toomre disks. 
We then adjusted the velocity scale to minimize the residuals between the data and our 
model rotation curves, finding (GM/R d ) 1 ' 2 = 584 & 540 km s" 1 for the maximum disk 
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Fig. 17. — The "maximum disk" and "control" simulations scaled to NGC 3198. The data 
for NGC 3198, shown in circles with error bars, are from van Albada et al. (1985). The 
solid lines show the full rotation curves, the dash-dotted lines the disk contribution and the 
dashed lines the halo contribution. The dotted line is the maximum disk as defined by van 
Albada et al. 



and control runs, respectively. The resulting scaled rotation curves are shown in Figure 
17, which also shows the maximum-disk fit of van Albada et al. As usual, both models fit 
NGC 3198 very well, which is another instance of the rotation curve degeneracy. Note that 
both systems are quite disk dominated, and that the maximum disk of van Albada et al. is 
perhaps even more disk dominated than our "maximum disk" simulation. 

Choosing length and velocity scales determines the time unit. With these adopted 
values, the duration of many of our experiments, 2 000 dynamical times, is equivalent to 
~ 10 Gyr. 



8. Conclusions 

8.1. Summary of Principal Results 

We have shown that the severe braking of the bar by dense isotropic halos reported in 
Paper I also occurs for other non-rotating, or backwards rotating, halos of similar density, 
whatever the shape of the halo velocity ellipsoid. In all such cases, the bar in the disk slows 
unacceptably in a few rotations. The bars in all our experiments with strongly prograde 
rotation in the halo were not braked as severely; the halo spins strongly in those cases for 
which friction ceased before 1Z rose above 1.4. 
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The existence of strong friction is in agreement with the theoretical prediction by 
Weinberg (1985), with earlier fully self-consistent simulations by Sellwood (1980), using 
a coarse grid, and by Athanassoula (1996) using a different iV-body method, and other 
work (e.g. Hernquist & Weinberg 1992). We also find that all our bars slow down as they 
lose angular momentum - a non-trivial result since bars are not rigid objects and could 
conceivably spin up (e.g. as a binary star) as angular momentum is lost. Remarkably, the 
bars in some of our experiments survived a truly drastic slow-down - more than a factor 
of five in many cases - providing further evidence that bars are in fact dynamically very 
robust objects (Miller Sz Smith 1979; Sellwood & Wilkinson 1993). We find no evidence to 
support the idea (Kormendy 1979) that strong braking of a bar might cause it to dissolve. 

In some of our simulations, the bar did a great deal of work on the halo - e.g. in the 
canonical run the halo gained 40% of the angular momentum of the disk. Nevertheless, the 
change in the halo density profile was quite minor (Figure |l|b). This example emphasizes 
that it is extremely difficult to change the density profile of a halo using interactions with 
the disk. 

Friction is generally reduced by lowering the density of the halo, and bars in maximum 
disks are able to remain fast (though only barely so) for large numbers of dynamical times, 
even in extensive, non-rotating halos, as reported in Paper I. The bar in our maximum disk 
simulation with an extensive halo is continuing to slow down even after 2000 dynamical 
times (Figure [Uja), suggesting that 71 might continue to rise. When scaled to NGC 3198, 
this continued evolution is too slow to matter, since we have followed it for 10 Gyr. But 
dynamical times are shorter in more luminous galaxies and our computed evolution lasts 
the equivalent of 5 Gyr in a galaxy with V max ~ 300 km s _1 . 

Friction is caused by a lag between the bar and an m = 2 distortion in the halo (Figure 
[5]). None of our halos is rotating sufficiently rapidly to be itself unstable to bar-forming 
modes (e.g. Sellwood & Valluri 1997), so the halo distortions are responses forced by the 
rotating bar in the disk, as is usually the case in dynamical friction. It is encouraging 
that we have been able to find some suggestion of the expected velocity dependence of the 
frictional force in our very crude analysis (Figure |15|), which suggests that the qualitative 
effect of halo rotation is predictable. 

Braking ceases once the forced response in the halo rotates in alignment with the bar in 
the disk. The gradual trapping of halo orbits into the driven non-axisymmetric potential is 
itself one of the principal sources of dynamical friction. It seems reasonable that trapping of 
halo orbits should involve less angular momentum loss for the bars in halos with an excess 
of particles with J z > 0, as we have found. Note that we have observed the locking of the 
halo response into alignment with the bar only in models with strongly confined halos. 
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8.2. Implications for barred galaxies 

As noted in the introduction, the rather small number of actual measurements of real 
barred galaxies all lie in the range 1 < 1Z < 1.4; the existence and locations of dust lanes 
in bars is indirect evidence that these low values pertain more generally. Thus our results 
require either that most strong bars are (1) too young to have been slowed significantly (2) 
exist in strongly rotating halos, or (3) are not embedded in dense halos. We review each of 
these in turn. 

If disks are not maximal, and halos do not rotate strongly, then bars must indeed be 
young objects to have remained fast today. Since the rate of bar slow-down scales with the 
halo density, the larger the required density, the younger they must be to avoid any slow 
cases. There is a suggestion that bars were rare in the early universe (Abraham et al. 1999), 
but they have certainly existed, probably in about their present numbers, since a redshift of 
one half. Their ages are therefore ^ 4 Gyr, or <: 800 dynamical times when our simulations 
are scaled to NGC 3198, which is plenty long enough for friction to have slowed the bars 
significantly, although perhaps not completely. 

Bars in galaxies which are significantly sub-maximal can remain fast for cosmologically 
interesting times only if the halo is anisotropic and rotates strongly throughout most of the 
disk region. The required halo angular momentum is very large, however. If all halo angular 
momentum arises from tidal torques in the early universe, the required value of A would 
be reached only rarely (Barnes & Efstathiou 1987; Steinmetz & Bartelmann 1995). Since 
^ 50% of all HSB disks containing strong bars (Eskridge et al. 2000), the vast majority 
cannot have A large enough to avoid bar slow down. 

Alternatively, one could imagine the inner halo to have been torqued up by some 
means. Tremaine & Ostriker (1999) suggest two ways to transfer angular momentum from 
the disk to the inner halo for precisely this purpose. We have found that if the halo has 
moderate central density, then it must have a high degree of rotation - fully half that of 
the disk, out to well beyond the bar's end. If the nearest strongly barred galaxy, our own 
Milky Way, has a sub-maximal disk, we require substantial rotation in the halo within the 
Solar radius for the bar to have avoided strong braking. Most torquing mechanisms would 
act equally on both the dark halo and any associated stars. Thus the absence of significant 
rotation in, for example, the metal-weak globular cluster system of the Milky Way (Harris 
2000) also suggests that the inner halo lacks the required angular momentum. 

We therefore conclude that bars in real galaxies remain fast because disks are maximal. 
Weiner et al. (2000) reach a similar conclusion on quite different grounds in the case of the 
barred galaxy NGC 4123. 
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8.3. General Implications 

It is often argued (e.g. Courteau & Rix 1999) that barred galaxies have lower density 
halos than do unbarred. This prejudice stems from the paper by Ostriker & Peebles (1973), 
who suggested that the only way to inhibit bar formation in a galaxy was to immerse the 
cool disk in a massive dynamically hot component. Not only is this argument fallacious 
(Toomre 1981; Sellwood & Evans 2000), but we here present evidence that bright barred 
galaxies have similar DM fractions as do their unbarred counterparts. 

In Paper I, we argued against the hypothesis that barred galaxies have less DM than 
those of the unbarred family: If the DM content varies continuously between maximum 
disk, fast bar, SB galaxies and massive halo SA galaxies, there should be many galaxies of 
intermediate dark matter content. Any strong bars that may form in such galaxies would 
therefore be fiercely braked, as in our experiments. Since no slow bars are known in HSB 
galaxies, we conclude that, either the distribution of dark matter is bimodal, or that all 
galaxies with moderate halo density have somehow avoided forming bars, both of which 
seem very unlikely, or that no galaxies are dark matter dominated. 

Tidal triggering can induce a bar in a galaxy that is stable when isolated (e.g. Noguchi 
1987). Such bars would be strongly braked if the target galaxy had been stabilized by a 
massive halo. The absence of known slow bars again argues against massive halos, but only 
weakly; if the rate of bar-inducing tidal interactions is low, then the measured sample may 
be simply too small to include a slow case. 

Empirical evidence against a systematic difference between barred and unbarred 
galaxies was presented by Bosma (1996) and more can be found in the data from 
Mathewson & Ford (1996). We use the apparent magnitudes in the I-band, recession 
velocities and V m given in their table, convert to absolute magnitudes assuming Hubble 
distances (for H Q = 60 km s _1 Mpc -1 ) and plot the Tully-Fisher relation for the 2368 
galaxies in their sample having recession speeds > 1, 000 km s _1 (to avoid absurdly 
inaccurate Hubble distances) in Figure [l8|(a). The line is fitted to all the data, but the 
2219 points are for "unbarred" galaxies and the crosses mark the 149 galaxies which 
Mathewson & Ford designate as barredfj The histograms in Figure |18|(b) show the 
distributions of velocity residuals about the fitted line for the barred and unbarred galaxies 
separately (scaled so that the area under each histogram is unity), suggesting a small 
offset in the sense that the barred galaxies have slightly lower V m at a given brightness. A 



4 Their sample excluded galaxies designated as barred but, as always happens, bars were identified after 
the data were taken. It is unclear whether these are typical bars, but since their barred fraction is extremely 
low, it seems likely that they flagged only the blatant, i.e. strong, bars. 
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Fig. 18. — The Tully-Fisher relation (a) and histograms of velocity residuals (b) for the 
sample collected by Mathewson & Ford (1996). Unbarred galaxies are plotted as points in 
(a), barred galaxies are marked by crosses and the slope and intercept of the fitted line in 
(a) are —0.126 & —0.566 respectively. The two histograms in (b) have been scaled so as to 
have equal (unit) area. 
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Kolmogorov-Smirnov test indicates that there is 3.5% probability that these two samples 
were drawn from the same parent distribution, suggesting a possibly significant difference. 
However, the offset disappears if we discard all galaxies fainter than Mi ~ —21, indicating 
that it arises from the faint galaxies only, as is apparent in Figure |T3(a). We conclude 
that there is no evidence here for a deficiency in DM content, relative to the unbarred 
galaxies, in the (few) bright barred galaxies in the Mathewson & Ford sample. Further 
Tully-Fisher work with properly constructed samples of barred/unbarred galaxies to confirm 
this conclusion would be highly desirable. 

We conclude that all bright HSB disk galaxies, barred or unbarred, are maximum disks. 
Supporting evidence for this conclusion is reviewed by Sellwood (1999). 

We also predict that any barred galaxy having a moderately dense halo should have a 
slow bar. Prime candidates to test this prediction may be found amongst galaxies believed 
to have significant DM fractions in their inner regions: the low luminosity galaxies (e.g. 
Persic & Salucci 1988; see Sellwood 1999) and low surface brightness galaxies (LSBs, 
Bothun et al. 1997; de Blok & McGaugh 1997). Bars in these systems are less common, but 
not unknown.^ If a strongly barred low-luminosity or LSB galaxy has even a moderately 
dense DM halo, it should have a high value of 1Z. Unfortunately, there are no reliable 
measurements of pattern speeds in such galaxies to test the prediction at this time. 
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critical read of the manuscript. The comments of the referee, Lia Athanassoula, were also 
helpful. This work was supported by NSF grant AST 96/17088 and NASA LTSA grant 
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A. Determination of an Equilibrium Distribution Function 

We here describe the creation of our equilibrium halo models in the presence of a 
massive disk. The procedure is identical to that followed by Raha et al. (1991), but has 
not previously been explained in detail in a published article. Because the equilibrium 
generated by this procedure is exact, there is no need to compute the evolution of the halo 
while it adjusts to an equilibrium from an approximate set-up (e.g. Barnes 1992; Hernquist 



We adopt the iterative approach to finding a DF first proposed by Prendergast & 
Tomer (1970) and developed by Jarvis & Freeman (1985) for two component systems 
and Kuijken & Dubinski (1995) for a three component system. Unlike these authors, 
however, we solve for the gravitational potential using the same numerical procedure 
that is used in the simulations, thereby incorporating any numerical idiosyncrasies of the 
potential determination into the solution for the DF; this strategy ensures that the particle 
distribution is in perfect equilibrium at the outset. 

We first choose a functional form for the DF of the halo 



where C is a normalization constant and T can be more or less any reasonable function of 
the classical isolating integrals, /„. In our axisymmetric potential, these are E and J z . The 
form of T adopted determines the density profile and shape of the resulting halo; functions 
of E alone tend to produce almost spherical halos (the disk makes them slightly oblate), 
while adding a dependence on J z generally produces more strongly spheroidal halos. 

A first approximation to the halo density p\ (R, z) can be determined from 



using any reasonable initial guess for the gravitational potential <&i(R,z). We assign mass 
to the grid to represent the smooth function pi(R, z), add the mass distribution of a smooth 
disk and solve for a new gravitational field $ 2 (R,z). We then determine p 2 (R,z) using 
the improved potential in (A2), and iterate until the potential distribution converges. The 
value of C can be adjusted at each iteration step to drive the solution to the desired halo 
mass. We find that the solution converges rapidly and that 10 iterations are usually ample. 

Note that although the halo density profile, and therefore the net rotation curve, 
cannot be specified in advance, the rapid convergence permits many models to be explored 
(for different mass ratio, truncation radius, choice of JF, etc.) from which one having the 
desired properties may be selected. 



1993). 



f = CF(h,I 2 ,...), 



(Al) 




(A2) 
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B. Quiet Start Procedures 



Since an TV-body simulation amounts to a numerical solution of the coupled collisionless 
Boltzmann and Poisson equations by the method of characteristics, it is clearly desirable 
to select the characteristics to be followed with care. Selecting particles at random from 
a DF (e.g. Hernquist, Spergel & Heyl 1993; Kuijken & Dubinski 1995) leads to \/jV-type 
variations in the number of particles generated in any given range of the integrals; in effect 
the model will have the dynamical properties of one with a DF slightly different from 
that intended, which has many disadvantages, especially when attempting comparisons 
with theoretical work. The following quiet start procedures lead to much higher quality 
simulations (and are also more efficient); every part has been described in some other 
publication, but for ease of reference we summarize them here. 

Strategies for the optimal selection of points are exactly analogous to those for 
the selection of abscissae in the numerical evaluation of multi-dimensional integrals. In 
that case, accuracy is improved whenever the dimensionality can be reduced by analytic 
integration over some of the coordinates. In our problem, we know the DF to be independent 
of orbital phases, since they must be uniformly populated in any equilibrium model. Note 
that, except when the DF is expressed in terms of actions, the density of particles in the 
sub-space of the integrals is not simply given by the DF; it needs to multiplied by a "density 
of states" function, which is the phase-space volume per unit interval of E and J z (Binney 
& Tremaine 1987, §4.4.5) 

For a razor-thin disk, the density of particles having energy E and angular momentum 
J z is (Sellwood & Athanassoula 1986) 



where / is the usual phase space density and t(E, J z ) is the period of one complete radial 
oscillation of a particle with the given (E, J z ). The latter generally has to be determined 
numerically. For a sphere with f(E,L) 



A/" disk (£, J z ) = 2nf(E, J z )t(E, J z ), 



(Bl) 



A/" sp here(-E, L) 



8n 2 Lf(E,L)r(E,L), 



(B2) 



(Binney & Tremaine 1987, problem 4-22) 



while, for a spheroid with f(E, J z ), 



jV S pheroid(£, Jz) = 47T 2 f(E, J Z )S(E, J z ) 



(B3) 



(Sellwood & Valluri 1997). In this last formula, S(E, J z ) is the cross-sectional area in (R, z) 
of the bounding torus in the meridional plane (Binney & Tremaine 1987, §3.2.1) and is 
easily evaluated numerically for arbitrary potentials. 
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We proceed by slicing accessible (E, J z ) space into a number, j = n E rij, of small areas 
in such a way that the integral of the DF over each area encloses a fraction 1/j of the total 
active mass; we generally choose n E ^> rij. We then select a point within each of these 
areas to determine the (E, J z ) values for an orbit. We avoid a regular raster of such points 
in (E, J z ) space while maintaining the desired smoothness, as follows: For every slice in 
energy, we choose rij equal spaced values of J z from the distribution of N\e, with the first 
value only determined as a random sub-fraction, and then select an E value within each 
area at random from the distribution of / M(E, J z )dJ z . 

Each (E, J z ) pair selected in this way specifies an orbit and we must next choose phases 
to determine both the initial position and velocity components of the particles. In contrast 
to the selection of integrals, experience suggests that the behavior of the model is much 
less sensitive to the manner in which some orbital phases are selected. In general, random 
selection from the appropriate distribution is adequate, although it is easy to improve upon 
random when desired for a particular application. Two examples are: Sellwood (1983) 
found it desirable to space several otherwise identical particles equally around a ring when 
searching for small- amplitude non-axisymmetric instabilities and Sellwood (1997) was able 
to quieten radial pulsations of a stable spherical model by spacing particles equally in radial 
phase. 

In a razor-thin disk, or in a sphere, the orbit lies in a plane in which particles oscillate 
between peri- and apocenter with full period t(E,L). The radial phases must be uniformly 
distributed, but the probability of selecting a particular radius varies inversely with the 
(non-uniform) radial speed. It is easiest to select a fraction of the radial period and 
integrate the orbit (usually numerically) for this time to determine the radius. The radial 
and azimuthal speeds are completely determined (except for the sign ambiguity of the 
radial speed) by the selected values of E, J z and r. The azimuthal phase and, for a sphere, 
the orientation of the plane can be selected in a straightforward manner. 

In spheroidal models, as here, the two classical integrals confine the particle to a torus 
in real space. When the desired DF does not depend upon a third integral, the probability 
density distribution for any given orbit is uniform in the meridional plane within the 
boundaries of the torus and selection of (R, z) pairs is straightforward. The values of E 
and J z almost determine the velocity components at the chosen position - all that remains 
is to direct the velocity component in the meridional plane, which should be uniformly 
distributed in angle. 
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C. Halo Angular Momentum 

In a general axisymmetric system, the dependence of the density p on J z is only 
through the even part of the DF (Lynden-Bell 1962). Thus angular momenta about the 
symmetry axis z can be reversed at random without affecting the equilibrium of the system. 
Changing the net angular momentum may, however, alter the stability of the system; in 
particular, Kalnajs (1977) has cautioned that a discontinuity in the DF at J z = can 
aggravate instabilities. We therefore adopt a scheme which ensures a smooth DF. 

If a prograde halo is desired, we define p a (x), where x = — J z j J z , max , to be the 
probablity of changing the sign of J z of a retrograde particle in a halo in which the 
maximum possible angular momentum at the truncation energy is J Z;max - (To generate a 
retrograde halo, we flip only particles having positive J z with probablity p a (x), where now 

X J z / Jz,max-) 

In order to make the DF continuous at J z = 0, we require p a — > as x — > 0. By 
extending the definition of p a (x) to be an odd function when x < 0, we can write the new 
distribution function as f(E, J z ) = f(E)[l — p a (x)]. 

We adopted the shifted, normalized Fermi function and its inverse, which both have 
the desired properties: 

| - [(e ax + I)- 1 - ±] / [(e- a + l)- 1 - \] large L 2 , h 

Pa(X) = { -J In ({* [(e- + I)- - \] + \Y l - l) small L z , h (C1) 

When a = 0, both functions are po(x) = x, which leads to a certain total L z ^. If the desired 
total L Zt h is greater (smaller) than this value, we use the large (small) expression and adjust 
a to yield the desired net angular momentum. 

In order to generate a radial variation in the net rotation, we make the probability a 
function of energy as follows: 

(ax — I — I)- 1 _ 1 

»<*> - + l r > A < C2 » 

with 

1 E < Ei = E min + (E max — E m [ n ) j 

N(E) = I 2s 3 - 3s 2 + 1 E > E > E 1 (C3) 

E > Eq = -Emin + (-Emax — -E'min)(^ + ^j) 

Here s = (E — Ei)/(E — -Ei), while i? m i n and -E max are constants of the system. We fix 
a = 30 to ensure the inner halo rotates strongly and vary d to adjust the energy at which 
strong rotation gives over to no rotation. 



